A Primer on Resurgent Transseries and Their Asymptotics -- 1802.10441 -- Physics Reports 809 (2019) 1-135

[Notebook by Roberto Vega (2022)]

Resurgent Analysis of a Quartic-Potential Integral

Phase Diagram from Stokes Phenomena

Figure 1

(*ContourPlot plots the gradient of the function. In the first case we chose 1 layer (bigger and smaller than 0),
while in the second case we chose 10 layers, to explain how it works*)
(*run OK*)
ContourPlot[Re[(1/2)(z^2)-(λ/24)(z^4)/.{λ->1,z->a+I b}],{a,-6,6},{b,-6,6},Contours->1,PlotLegends->Automatic]
ContourPlot[Re[(1/2)(z^2)-(λ/24)(z^4)/.{λ->1,z->a+I b}],{a,-6,6},{b,-6,6},Contours->10,PlotLegends->Automatic]

Figure 2

(*run OK*)
V[z_,x_]:=(1/2)(z^2)-(x/24)(z^4);
(*The module function defines internal names such that we don't have issues with defining variables*)
(*The next piece of code computes the minima of the potential in terms of x*)
minima[x_]:=Module[{z},
z/.Solve[D[V[z,x],z]==0,z]
]
minima[x]

(*runs SLOW*)
Module[{i,θl,x},
θl={0,(π/16),((7π)/16),(π/2),((9π)/16),((15π)/16),π,((17π)/16),((23π)/16),((3π)/2),((25π)/16),((31π)/16)};
For[i=1,i≤Length[θl],i++,
(*For each of the angles we are interested in checking, we define x to be the exponential of this angle*)
x=Exp[I θl[[i]]];
(*The plot we will generate depends on which saddles are dominating the integral for the specific x at hand.
The only thing that will change between the three cases in the Which function are which ContourStyles are dashed (for subleading saddles) and which ones are solid (for leading saddles). Also, we will plot the saddle themselves using Circle[] for the leading and Disk[] for the subleading.*)
Which[
(*This first case corresponds to the two symmetric minima being leading*)
Abs[Exp[V[minima[x][[1]],x]]]<Abs[Exp[V[minima[x][[2]],x]]],
Print[Show[
RegionPlot[Re[V[z,x]/.z->a+I b]>0,{a,-5,5},{b,-5,5},PlotStyle->LightOrange,PlotLabel->θl[[i]]],RegionPlot[Re[V[z,x]/.z->a+I b]<0,{a,-5,5},{b,-5,5},PlotStyle->LightBlue],ContourPlot[Im[(V[z,x]/.z->a+I b)-V[minima[x][[2]],x]]==0,{a,-5,5},{b,-5,5},ContourStyle->Directive[Red,Dashed]],ContourPlot[Im[(V[z,x]/.z->a+I b)-V[minima[x][[1]],x]]==0,{a,-5,5},{b,-5,5},ContourStyle->Black],ListPlot[{ReIm[minima[x][[2]]],ReIm[minima[x][[3]]]},PlotStyle->Red,PlotMarkers->{Graphics[{Red,Circle[]}],0.03}],ListPlot[{ReIm[minima[x][[1]]]},PlotMarkers->{Graphics[{Black,Disk[]}],0.03}]]],
(*This second case corresponds to the two 0 minimum being leading*)
Abs[Exp[V[minima[x][[1]],x]]]>Abs[Exp[V[minima[x][[2]],x]]],
Print[Show[
RegionPlot[Re[V[z,x]/.z->a+I b]>0,{a,-5,5},{b,-5,5},PlotStyle->LightOrange,PlotLabel->θl[[i]]],RegionPlot[Re[V[z,x]/.z->a+I b]<0,{a,-5,5},{b,-5,5},PlotStyle->LightBlue],ContourPlot[Im[(V[z,x]/.z->a+I b)-V[minima[x][[2]],x]]==0,{a,-5,5},{b,-5,5},ContourStyle->Red],ContourPlot[Im[(V[z,x]/.z->a+I b)-V[minima[x][[1]],x]]==0,{a,-5,5},{b,-5,5},ContourStyle->Directive[Black,Dashed]],ListPlot[{ReIm[minima[x][[2]]],ReIm[minima[x][[3]]]},PlotStyle->Red,PlotMarkers->{Graphics[{Red,Disk[]}],0.03}],ListPlot[{ReIm[minima[x][[1]]]},PlotMarkers->{Graphics[{Black,Circle[]}],0.03}]]],
(*This third case corresponds to the two all minima having the same weight*)
Abs[Exp[V[minima[x][[1]],x]]]==Abs[Exp[V[minima[x][[2]],x]]],
Print[Show[
RegionPlot[Re[V[z,x]/.z->a+I b]>0,{a,-5,5},{b,-5,5},PlotStyle->LightOrange,PlotLabel->θl[[i]]],RegionPlot[Re[V[z,x]/.z->a+I b]<0,{a,-5,5},{b,-5,5},PlotStyle->LightBlue],ContourPlot[Im[(V[z,x]/.z->a+I b)-V[minima[x][[2]],x]]==0,{a,-5,5},{b,-5,5},ContourStyle->Red],ContourPlot[Im[(V[z,x]/.z->a+I b)-V[minima[x][[1]],x]]==0,{a,-5,5},{b,-5,5},ContourStyle->Directive[Black]],ListPlot[{ReIm[minima[x][[2]]],ReIm[minima[x][[3]]]},PlotStyle->Red,PlotMarkers->{Graphics[{Red,Disk[]}],0.03}],ListPlot[{ReIm[minima[x][[1]]]},PlotMarkers->{Graphics[{Black,Disk[]}],0.03}]]]
]
]
]

Basic Formulae for Resurgent Analysis

(*run OK*)
(*In here, we have already set hbar=2π*)
Z0[x_]:=(Sqrt[((2π)/(2π))])Sum[((((2/3)))^n)  ((Factorial[4n])/(2^(6n)*Factorial[2n]*Factorial[n]))(x^n),{n,0,120}]

(*For readability and in order to see the Borel transform explictly, we define the coefficients of Z0*)
coefZ0[n_]:=(Sqrt[((2π)/(2π))])((((2/3)))^n)  ((Factorial[4n])/(2^(6n)*Factorial[2n]*Factorial[n]))
(*The Borel transform is thus given by*)
BorelZ0[s_]:=Sum[((coefZ0[n])/(Gamma[n+1]))(s^n),{n,0,120}]

Z1[x_]:=((-I)/(Sqrt[2]))(Sqrt[((2π)/(2π))])Sum[(((-1))^n)((((2/3)))^n)  ((Factorial[4n])/(2^(6n)*Factorial[2n]*Factorial[n]))(x^n),{n,0,120}]
coefZ1[n_]:=((-I)/(Sqrt[2]))(Sqrt[((2π)/(2π))])(((-1))^n)((((2/3)))^n)  ((Factorial[4n])/(2^(6n)*Factorial[2n]*Factorial[n]))
BorelZ1[s_]:=Sum[((coefZ1[n])/(Gamma[n+1]))(s^n),{n,0,120}]

(*The Padé Approximant has been taken for the variable s, around 0, and taking 60 as the order of the polynomial both in the numerator and the denominator. 
If one were to play with different order polynomials in the numerator and denominator one would write, e.g.,
PadeApproximant[BorelZ0[s],{s,0,{numdeg,dendeg}}]
Also, the fact that we write PadeZ0[s_]:=PadeZ0[s]= is because we want Mathematica to save the result of each computation, making future calculations much faster*)

PadeZ0[s_]:=PadeZ0[s]=PadeApproximant[BorelZ0[s],{s,0,60}]
PadeZ1[s_]:=PadeZ1[s]=PadeApproximant[BorelZ1[s],{s,0,60}]
(*Example of chosing different degrees*)
PadeApproximant[BorelZ0[s],{s,0,{3,6}}]

(*Run OK*)
(*Now, the resumed Z function is obtained by integrating the Padé approximated Borel transforms.
For computational efficiency it is better to compute the Padé at a fixed variable s (which will be saved into Mathematica's memory) and substitute the value of s in the memorised expression rather than computing this expression for each value*)
Z[σ0_,σ1_,absx_,θ_]:=Z[σ0,σ1,absx,θ]=σ0*NIntegrate[(E^((-s)/( absx Exp[I*θ])))*PadeZ0[s]/.s->auxs,{auxs,0,Exp[I*θ]*Infinity}]+σ1*Exp[((-3)/(2 absx Exp[I*θ]))]*NIntegrate[(E^(((-s)/( absx Exp[I*θ])) ))*PadeZ1[s]/.s->auxs,{auxs,0,Exp[I*θ]*Infinity}]
(*As a check we compute the residues of the Padé approximation, i.e., the zeroes of the denominator to see where the branch cut of the function lies on the complex Borel plane.*)
Module[{polesofPadeZ0,polesofPadeZ1},
polesofPadeZ0=Solve[Denominator[PadeZ0[s]]==0];
Print[Show[
ListPlot[Table[ReIm[s/.t],{t,polesofPadeZ0}],PlotRange->{{0,10},Automatic},AspectRatio->1/3],
ListPlot[{ReIm[(3/2)]},PlotStyle->Red],ImageSize->Large]];
polesofPadeZ1=Solve[Denominator[PadeZ1[s]]==0];
Print[Show[
ListPlot[Table[ReIm[s/.t],{t,polesofPadeZ1}],PlotRange->{{-10,0},Automatic},AspectRatio->1/3],
ListPlot[{ReIm[-(3/2)]},PlotStyle->Red],ImageSize->Large]];

]

(*DO NOT RUN in real time*)

(*In here, we want to see the monodromy of the function Z that we have resummed. We have implemented by hand the Stokes transitions for the σ0, σ1 transseries parameters, while we have varied the angle of x, denoted by θ, from 0 to 8π*)
Print["The order of the colors is Cyan->Blue->Purple->Red->Orange->Yellow->Green->Pink"]
Show[
ListPlot[Table[ReIm[Z[1,-2,0.21,θ]],{θ,0,π,0.01π}],PlotStyle->Cyan,PlotRange->All,ImageSize->Full],ListPlot[Table[ReIm[Z[-1,-2,0.21,θ]],{θ,π,2π,0.01π}],PlotStyle->Blue,PlotRange->All],
ListPlot[Table[ReIm[Z[-1,0,0.21,θ]],{θ,2π,3π,0.01π}],PlotStyle->Purple,PlotRange->All],
ListPlot[Table[ReIm[Z[-1,0,0.21,θ]],{θ,3π,4π,0.01π}],PlotStyle->Red,PlotRange->All],
ListPlot[Table[ReIm[Z[-1,2,0.21,θ]],{θ,4π,5π,0.01π}],PlotStyle->Orange,PlotRange->All],
ListPlot[Table[ReIm[Z[1,2,0.21,θ]],{θ,5π,6π,0.01π}],PlotStyle->Yellow,PlotRange->All],
ListPlot[Table[ReIm[Z[1,0,0.21,θ]],{θ,6π,7π,0.01π}],PlotStyle->Green,PlotRange->All],
ListPlot[Table[ReIm[Z[1,0,0.21,θ]],{θ,7π,8π,0.01π}],PlotStyle->Pink,PlotRange->All],
PlotRange->All
]
(*We have also restricted the plot to see the continuity of the smaller results*)
Print["The order of the colors is Cyan->Blue->Purple->Red->Orange->Yellow->Green->Pink"]
Show[
ListPlot[Table[ReIm[Z[1,-2,0.21,θ]],{θ,0,π,0.01π}],PlotStyle->Cyan,PlotRange->All,ImageSize->Full],ListPlot[Table[ReIm[Z[-1,-2,0.21,θ]],{θ,π,2π,0.01π}],PlotStyle->Blue,PlotRange->All],
ListPlot[Table[ReIm[Z[-1,0,0.21,θ]],{θ,2π,3π,0.01π}],PlotStyle->Purple,PlotRange->All],
ListPlot[Table[ReIm[Z[-1,0,0.21,θ]],{θ,3π,4π,0.01π}],PlotStyle->Red,PlotRange->All],
ListPlot[Table[ReIm[Z[-1,2,0.21,θ]],{θ,4π,5π,0.01π}],PlotStyle->Orange,PlotRange->All],
ListPlot[Table[ReIm[Z[1,2,0.21,θ]],{θ,5π,6π,0.01π}],PlotStyle->Yellow,PlotRange->All],
ListPlot[Table[ReIm[Z[1,0,0.21,θ]],{θ,6π,7π,0.01π}],PlotStyle->Green,PlotRange->All],
ListPlot[Table[ReIm[Z[1,0,0.21,θ]],{θ,7π,8π,0.01π}],PlotStyle->Pink,PlotRange->All],
PlotRange->{{-0.25,0.25},{-0.25,0.25}}
]

Asymptotics and Large-Order Behaviour: Partition Function

Figure 6

(*Run OK*)
coeffZ0[n_]:=(1/((Sqrt[2])π))((((2/3)))^n) Sum[(((-1))^k) (((4k)!(n-k-1)!)/((2^(6k))k!(2k)!)),{k,0,n-1}]
(*The next function is the n-th Richardson transform of the series s*)
rt0[n_,s_]:=rt0[n,s]=If[n≥Length[s],Return[{}],
Module[{out=Array[0,Length[s]-n]},
(If[n==0,out=s,Module[{k=1},
For[k=1,k≤Length[s]-n,k++,
out[[k]]=rt0[n-1,s][[k+1]]+(k/n)(rt0[n-1,s][[k+1]]-rt0[n-1,s][[k]])]]];Return[out])]]
(*Since even and odd coefficients tend to the same quantity at different paces, Richardson transforms won't work unless we consider only the even or only the odds.*)
t=Table[((coeffZ0[2n]2π(3/2)^(2n))/(Gamma[2n])),{n,0,100}];
richardsont=rt0[5,t];

Show[
ListPlot[t,PlotStyle->Red],
ListPlot[richardsont,PlotStyle->Blue],
PlotRange->{Automatic,{1.390,1.420}}]

First Steps Towards Nonlinear Resurgent Analysis

Figure 7

(*Run OK*)
(*This is one way of defining coefficients recursively. Also, the most elegant way, we will present later another possibility.*)
F0[n_]:=F0[n]=If[n==1,(1/8),
(2/3)((n-1) F0[n-1]+(1/n)(∑_(k=1)^(n-2))k(n-1-k)F0[k]F0[n-1-k])
]
(*Φ0 is the Borel transform of F0*)
Φ0[s_]:=(Sqrt[(1/(2π))])Sum[((F0[n])/(Gamma[n]))(s^(n-1)),{n,1,121}]
PadeΦ0[s_]:=PadeΦ0[s]=PadeApproximant[Φ0[s],{s,0,60}]

polesPadeΦ0=s/.Solve[Denominator[PadeΦ0[s]]==0];

Show[
ListPlot[ReIm[polesPadeΦ0],PlotStyle->PointSize[0.03]],
ListPlot[{{(3/2),0}},PlotStyle->{Black, PointSize[0.03]}],
PlotRange->{{0,6},{-0.6,0.6}}]

(*Run OK*)
F1[n_]:=F1[n]=Which[
n==1,-(1/4),
n==2,(1/32),
True,
-(2/(3(n)))((n-1)*n*F1[n-1]+2Sum[(n-1-l)*l*F1[n-1-l]F0[l],{l,1,n-2}]+3*Sum[(n-l)*F0[n-l]F1[l],{l,1,n-1}])
]

Φ1[s_]:=Sum[((F1[n])/(Factorial[n]))*s^n,{n,1,120}]
PadeΦ1[s_]:=PadeΦ1[s]=PadeApproximant[Φ1[s],{s,0,60}]

polesPadeΦ1=s/.Solve[Denominator[PadeΦ1[s]]==0];

Show[
ListPlot[ReIm[polesPadeΦ1],PlotStyle->PointSize[0.03]],
ListPlot[{{(3/2),0},{-(3/2),0}},PlotStyle->{Black, PointSize[0.03]}],
PlotRange->{{-6,6},{-0.6,0.6}}]

Asymptotics and Large-Order Behaviour : Free Energy

Figure 12

(*Run OK*)
rt0[n_,s_]:=rt0[n,s]=If[n≥Length[s],Return[{}],
Module[{out=Array[0,Length[s]-n]},
(If[n==0,out=s,Module[{k=1},
For[k=1,k≤Length[s]-n,k++,
out[[k]]=rt0[n-1,s][[k+1]]+(k/n)(rt0[n-1,s][[k+1]]-rt0[n-1,s][[k]])]]];Return[out])]]
rtrec[r_,n_,s_]:=rtrec[r,n,s]=If[n r+r≥Length[s],Return[{}],
Module[{out=Array[0,Length[s]-n r-r]},
(If[r==0,out=rt0[n,s],
Module[{newsuc=Array[0,Length[s]-n r-r+1]},
(For[k=1,k≤Length[s]-n r-r+1,k++,
newsuc[[k]]=(s[[k+1]])];
out=rt0[n,newsuc])]];Return[out])]]

A=3/2;
(*An alternative way of defining recursive coefficients is to accumulate them in a list as follows*)
F0in={(1/8),(1/12)};
Module[{k},
For[k=2,k≤124,k++,
(*The next piece of code appends the result of the second entry computation at the end of the list F0in*)
AppendTo[F0in,(2/3)(k*F0in[[k]]+(1/(k+1))Sum[(k-l)*l*F0in[[k-l]]F0in[[l]],{l,1,k-1}])]]
]

F010=((-I)/(Sqrt[2]));
F1in={-(1/4)F010,(1/32)F010};

Module[{k},
For[k=2,k≤124,k++,AppendTo[F1in,-(2/(3(k+1)))(k*(k+1)*F1in[[k]]+2Sum[(k-l)*l*F1in[[k-l]]F0in[[l]],{l,1,k-1}]+3*Sum[(k-l+1)*F0in[[k-l+1]]*F1in[[l]],{l,1,k}])]]
]


WeightedF0in=Table[-F0in[[i]]*(((2π *A^i)/(Factorial[i-1]))),{i,1,120}];
F0plot=Table[{i,WeightedF0in[[i]]},{i,1,120}];
RT1F0plot=Table[{i,rt0[1,WeightedF0in][[i]]},{i,1,Length[rt0[1,WeightedF0in]]}];
RT5F0plot=Table[{i,rt0[5,WeightedF0in][[i]]},{i,1,Length[rt0[5,WeightedF0in]]}];

Show[
ListPlot[{F0plot,RT1F0plot,RT5F0plot},PlotStyle->{Red,Green,Blue}],
Plot[-(Sqrt[2]),{x,0,120},PlotStyle->Purple],
PlotRange->{{0,120},{-1.420,-1.390}},
ImageSize->Large
]

(*Run OK*)
WeightedF0in1=Table[-F0in[[i]]*(((2π *A^i*i)/(Factorial[i-1])))+(Sqrt[2])i,{i,1,120}];
F0plot1=Table[{i,WeightedF0in1[[i]]},{i,1,120}];
RT1F0plot1=Table[{i,rt0[1,WeightedF0in1][[i]]},{i,1,Length[rt0[1,WeightedF0in1]]}];
RT5F0plot1=Table[{i,rt0[5,WeightedF0in1][[i]]},{i,1,Length[rt0[5,WeightedF0in1]]}];

Show[
ListPlot[{F0plot1,RT1F0plot1,RT5F0plot1},PlotStyle->{Red,Green,Blue}],
Plot[(3/8)(Sqrt[2]),{x,0,120},PlotStyle->Purple],
PlotRange->{{0,120},{0.52,0.56}},
ImageSize->Large
]

Check of equations (2.160), (2.166)

Chop[((N[-RT5F0plot[[115]][[2]]]I-I(Sqrt[2]))/(I(Sqrt[2])))]
Chop[((N[-RT5F0plot1[[115]][[2]]]I+((3I)/8)(Sqrt[2]))/(-((3I)/8)(Sqrt[2])))]

Figure 13: Implementing resummations

(*DO NOT RUN in real time*)
(*The ratio of Gamma functions can be obtained in Mathematica by doing Series expansions. However, it is very computationally inconvenient. We will use the expansion of the ratio of Gamma functions that NIST provides, which makes use of generalized Bell polynomials. These are encoded in Mathematica under the NorlundB function*)
G[k_,a_,b_]:=Binomial[a-b,k]NorlundB[k,a-b+1,a]
ratioGamma[k_,a_,b_,h_]:=(k^(a-b))Sum[((G[l,a,b])/(k^l)),{l,0,120-h}]
χ[0,1][k_]:=χ[0,1][k]=(1/(((1-0))^k))Sum[ratioGamma[k,-h,0,h]F1in[[h+1]]((((1-0)A))^h),{h,0,120}]

(*DO NOT RUN in real time*)
(*We take the first 100 coefficients since we have to do a Borel-Padé approximation of the series*)
coeffχ[0,1][g_]:=CoefficientList[χ[0,1][k]/.k->(1/z),z][[g]]
(*This is just a check on the growth of the coefficients. We sometimes use this to determine the specific growth of our coefficients, i.e., if we have to substract g!, (g+1)!, (2g)!, etc.*)

ListPlot[Table[((coeffχ[0,1][g+1])/(coeffχ[0,1][g](g+1))),{g,1,100}]]
RSolve[((f[g+1])/(f[g](g+1)))==A,f[g],g]//FullSimplify

(*DO NOT RUN in real time*)
Borelχ[0,1]=Sum[((coeffχ[0,1][g+1])/(Gamma[g+1]))(s^g),{g,0,120}];
Padeχ[0,1]=PadeApproximant[Borelχ[0,1],{s,0,60}];

Module[{poles},
poles=s/.Solve[Denominator[Padeχ[0,1]]==0];
ListPlot[Table[ReIm[t],{t,poles}]]
]
Resummedχ[0,1][k_]:=Resummedχ[0,1][k]=NIntegrate[SetPrecision[(E^(-t ))Padeχ[0,1]/.s->((t/k)),2000],{t,0,Exp[-(π/4)I]Infinity},MaxRecursion->50,WorkingPrecision->100]

leftfigure13list=Table[-(2^k)Re[Resummedχ[0,1][k]],{k,1,100}];

(*DO NOT RUN in real time*)
Show[
ListPlot[{8leftfigure13list,8rt0[5,leftfigure13list]},PlotStyle->{Red,Blue}],
Plot[1,{g,0,100},PlotStyle->Purple]
]

(*DO NOT RUN in real time*)

rightfigure13list=Table[SetPrecision[((Im[F0in[[k]]*(((2π I*(3A)^k)/(Gamma[k])))])/(-(3^k)(-8)*Im[Resummedχ[0,1][k]])),100],{k,1,120}]//N
Show[
ListPlot[{rightfigure13list,rt0[1,rightfigure13list]},PlotStyle->{Red,Green}],
Plot[1,{g,0,100},PlotStyle->Purple]
]

Exercise 89

u00=Table[0,{i,1,200}];
u00[[1]]=(I/(Sqrt[2]));
u00[[2]]=((1+(Sqrt[2]))/2)(1-I);
(*N[(I/(Sqrt[2]))((k+1)*(k+2)*u00[[k]]-(1+I(Sqrt[2]))*(k+2)*u00[[k+1]]+Sum[(2+I(Sqrt[2]))*u00[[h]]*u00[[k-h+1]]-I(Sqrt[2])(k-h+1)u00[[h]]u00[[k-h]],{h,1,k-1}]),1000]*)

For[k=1,k≤Length[u00]-2,k++,{u00[[k+2]]=N[(I/(Sqrt[2]))((k+1)*(k+2)*u00[[k]]-(1+I(Sqrt[2]))*(k+2)*u00[[k+1]]+Sum[(1+I(Sqrt[2]))*u00[[h]]*u00[[k-h+1]],{h,1,k}]-Sum[I(Sqrt[2])(k-h+1)u00[[h]]u00[[k-h]],{h,1,k-1}]+Sum[Sum[u00[[h]]u00[[f]]u00[[k-h-f]],{f,1,k-h-1}],{h,1,k-1}]),1000]
}]
Φ00[s_]:=Sum[((u00[[n]])/(Factorial[n]))*s^{n},{n,1,200}]
pade00[x_]:=N[PadeApproximant[Φ00[x],{x,0,100}],1000]
poles00=TransferFunctionPoles[TransferFunctionModel[pade00[x],x]][[1]][[1]];
poles00print=Table[0,{i,1,Length[poles00]}];

For[i=1,i≤Length[poles00],i++,poles00print[[i]]={Re[poles00][[i]],Im[poles00[[i]]]}];

Print[Show[ListPlot[poles00print,PlotRange->{{-6,6},{-0.6,0.6}},PlotStyle->PointSize[0.03]],ListPlot[{{1.5,0},{-1.5,0}},PlotRange->{{0,6},{-0.6,0.6}},PlotStyle->{Black, PointSize[0.03]}]]]

pade00[x_]:=N[PadeApproximant[Φ00[x][[1]],{x,0,100}],1000]
poles00=TransferFunctionPoles[TransferFunctionModel[pade00[x],x]][[1]][[1]];
poles00print=Table[0,{i,1,Length[poles00]}];

For[i=1,i≤Length[poles00],i++,poles00print[[i]]={Re[poles00][[i]],Im[poles00[[i]]]}];

Dots1=Table[{i+1,0},{i,1,6}];
Dots2=Table[{0,(i+1)(Sqrt[2])},{i,1,6}];
Dots3=Table[{i,j(Sqrt[2])},{i,1,6},{j,1,6}];
Print[Show[ListPlot[poles00print,PlotRange->{{-1,6},{-1,6}},ImageSize->1000,PlotStyle->PointSize[0.01]],ListPlot[{{1,0},{0,(Sqrt[2])}},PlotRange->{{0,6},{-0.6,0.6}},PlotStyle->{Green, PointSize[0.01]}],ListPlot[{Dots1,Dots2},PlotRange->{{0,6},{-0.6,0.6}},PlotStyle->{Orange, PointSize[0.01]}],ListPlot[Dots3,PlotRange->{{-1,6},{-0.6,0.6}},PlotStyle->Yellow,PlotMarkers->{Graphics[{Gray,Circle[]}],0.018}],ListPlot[{{0,0}},PlotRange->{{-1,6},{-0.6,0.6}},PlotStyle->Red]]]

We have to construct all the series of (m | n) with m <= 3 and n <= 2. Those are (0|0), (0|1),  (0|2), (1|0), (2|0), (3|0), (1|1), (1|2), (2|1), (2|2), (3|1)

num=50;
prec=1000;

Clear[u]
u[m_,n_,g_]:=u[m,n,g]=If[m<0||n<0||g<0,Return[0],
If[m==0&&n==0,
If[g==0,Return[N[(I/(Sqrt[2])),prec]],
If[g==1,Return[N[1+(I/(2(Sqrt[2]))),prec]],
Return[N[(1/(-I(Sqrt[2])))(u[0,0,g-1](-(g-1)(1+I(Sqrt[2])))+u[0,0,g-2]*((g-1)(g))+(1+I(Sqrt[2]))Sum[u[0,0,h]u[0,0,g-h-1],{h,0,g-1}]-Sum[u[0,0,h]u[0,0,g-h-2](I(Sqrt[2])(g-h-1)),{h,0,g-2}]+Sum[u[0,0,h]u[0,0,f]u[0,0,g-h-f-2],{h,0,g-2},{f,0,g-h-2}]),prec]]]],
If[m+n==1,
If [g==0,If[m==1,Return[N[1,prec]],Return[N[0,prec]]],
Return[N[(u[0,0,g-1](g)(g+1)+Sum[u[a,b,h]u[m-a,n-b,g-h-1]((1+I(Sqrt[2]))-I(Sqrt[2])((m-a)+(n-b)I(Sqrt[2]))),{h,1,g-1},{a,0,m},{b,0,n}]+Sum[u[a,b,0]u[m-a,n-b,g]((1+I(Sqrt[2]))-I(Sqrt[2])((m-a)+(n-b)I(Sqrt[2]))),{a,1,m},{b,1,n}]+Sum[u[a,b,g]u[m-a,n-b,0]((1+I(Sqrt[2]))-I(Sqrt[2])((m-a)+(n-b)I(Sqrt[2]))),{a,0,m-1},{b,0,n-1}]-Sum[u[a,b,h]u[m-a,n-b,g-h-1](I(Sqrt[2])(g-h)),{h,0,g-1},{a,0,m},{b,0,n}]+Sum[u[a,b,h]u[d,e,f]u[m-a-d,n-b-e,g-1-h-f],{h,0,g-1},{a,0,m},{b,0,n},{f,0,g-h-1},{d,0,m-a},{e,0,n-b}])/((1+I(Sqrt[2]))(g+1)-2(m+I(Sqrt[2])n)(g+1)+u[0,0,0](I(Sqrt[2])(m+I(Sqrt[2])n)-2(1+I(Sqrt[2])))),prec]]],
If[g==0,Return[N[0,prec]],
If[g==1,
Return[N[(u[m,n,0](2(m+I(Sqrt[2])n)-(1+I(Sqrt[2])))+Sum[u[a,b,0]u[m-a,n-b,0]((1+I(Sqrt[2]))-I(Sqrt[2])((m-a)+I(Sqrt[2])(n-b))),{a,0,m},{b,0,n}])/((1+I(Sqrt[2]))(m+I(Sqrt[2])n)-(m+I(Sqrt[2])n)^2-I(Sqrt[2])),prec]],
Return[N[(u[m,n,g-1](2(m+I(Sqrt[2])n)(g-1)-(1+I(Sqrt[2]))(g-1))+u[m,n,g-2]((g-1)(g))+Sum[u[a,b,h]u[m-a,n-b,g-h-1]((1+I(Sqrt[2]))-I(Sqrt[2])((m-a)+I(Sqrt[2])(n-b))),{h,0,g-1},{b,0,n},{a,0,m}]-Sum[u[a,b,h]u[m-a,n-b,g-h-2](I(Sqrt[2])(g-h-1)),{h,0,g-2},{b,0,n},{a,0,m}]+Sum[u[a,b,h]u[d,e,f]u[m-a-d,n-b-e,g-2-h-f],{h,0,g-2},{a,0,m},{b,0,n},{f,0,g-h-2},{d,0,m-a},{e,0,n-b}])/((1+I(Sqrt[2]))(m+I(Sqrt[2])n)-(m+I(Sqrt[2])n)^2-I(Sqrt[2])),prec]]]]]]]

Clear[Φ32]
Clear[pade32]
Φ32[s_]:=Φ32[s]=N[Sum[((u[3,2,g])/(Factorial[g]))*s^{g},{g,0,num}],prec]
pade32[x_]:=pade32[x]=N[PadeApproximant[Φ32[x][[1]],{x,0,(num/2)}],prec]
poles32=TransferFunctionPoles[TransferFunctionModel[pade32[x],x]][[1]][[1]];
poles32print=Table[0,{i,1,Length[poles32]}];
For[i=1,i≤Length[poles32],i++,poles32print[[i]]={Re[poles32][[i]],Im[poles32[[i]]]}];
Length[poles32print]
Φ32[s];

GDots=Delete[Table[{(-i+2),((-j+2)(Sqrt[2]))},{i,1,5},{j,1,5}],{{1,1},{5,4}}];
Dots1=Table[{i+1,0},{i,1,6}];
Dots2=Table[{i+1,-2(Sqrt[2])},{i,1,1}];
Dots3=Table[{0,(i+1)(Sqrt[2])},{i,1,6}];
Dots4=Table[{-2,(i+1)(Sqrt[2])},{i,1,1}];
Dots5=Table[{(-i+4),(-j+4)(Sqrt[2])},{i,1,8},{j,1,8}];
Print[Show[ListPlot[poles32print,PlotRange->{{-3.5,3.5},{-3.5,3.5}},ImageSize->1000,PlotStyle->PointSize[0.01]],ListPlot[GDots,PlotRange->{{0,6},{-0.6,0.6}},PlotStyle->{{Green, PointSize[0.01]},{Green, PointSize[0.01]}}],ListPlot[{Dots1,Dots2,Dots4},PlotRange->{{0,6},{-0.6,0.6}},PlotStyle->{{Orange, PointSize[0.01]},{Orange, PointSize[0.01]},{Orange, PointSize[0.01]}}],ListPlot[Dots5,PlotRange->{{-1,6},{-0.6,0.6}},PlotStyle->Yellow,PlotMarkers->{Graphics[{Gray,Circle[]}],0.018}],ListPlot[{{0,0}},PlotRange->{{-1,6},{-0.6,0.6}},PlotStyle->Red]]]

Φ11[s_]:=Sum[((u[1,1,n])/(Factorial[n]))*s^{n},{n,0,num}]
Φ11[s];
pade11[x_]:=N[PadeApproximant[Φ11[x][[1]],{x,0,(num/2)}],1000]
poles11=TransferFunctionPoles[TransferFunctionModel[pade11[x],x]][[1]][[1]];


poles11print=Table[0,{i,1,Length[poles11]}];
For[i=1,i≤Length[poles11],i++,poles11print[[i]]={Re[poles11][[i]],Im[poles11[[i]]]}];
Length[poles11print]

GDots=Delete[Table[{(-i+2),((-j+2)(Sqrt[2]))},{i,1,5},{j,1,5}],{{1,1},{5,4}}];
Dots1=Table[{i+1,0},{i,1,6}];
Dots2=Table[{i+1,-2(Sqrt[2])},{i,1,1}];
Dots3=Table[{0,(i+1)(Sqrt[2])},{i,1,6}];
Dots4=Table[{-2,(i+1)(Sqrt[2])},{i,1,1}];
Dots5=Table[{(-i+4),(-j+4)(Sqrt[2])},{i,1,8},{j,1,8}];
Print[Show[ListPlot[poles11print,PlotRange->{{-3.5,3.5},{-3.5,3.5}},ImageSize->1000,PlotStyle->PointSize[0.01]],ListPlot[GDots,PlotRange->{{0,6},{-0.6,0.6}},PlotStyle->{{Green, PointSize[0.01]},{Green, PointSize[0.01]}}],ListPlot[{Dots1,Dots2,Dots4},PlotRange->{{0,6},{-0.6,0.6}},PlotStyle->{{Orange, PointSize[0.01]},{Orange, PointSize[0.01]},{Orange, PointSize[0.01]}}],ListPlot[Dots5,PlotRange->{{-1,6},{-0.6,0.6}},PlotStyle->Yellow,PlotMarkers->{Graphics[{Gray,Circle[]}],0.018}],ListPlot[{{0,0}},PlotRange->{{-1,6},{-0.6,0.6}},PlotStyle->Red]]]